function [fval, eqm] = eqm_ss_iid_case_6(x, par, fun, eqm_all, cali_step)

% -----> Input:
% x:   initial v [L, B, ALPHA_0]
% par: parameters
% fun: functions to be used
% eqm_all: whether to compute all equilibrium outcomes; 0 or 1 values
% cali_step: 1 means general calibration; 2 means calibrating the entry cost

% -----> Output
% fval: residuals of equilibrium conditions
% eqm: equilibrium variables


L       = x(1); 
B       = x(2);
ALPHA_0 = x(3);

%%%% Controling the integration
int_method = par.int_method;
Rel_Tol    = par.Rel_Tol;
Abs_Tol    = par.Abs_Tol;

if cali_step == 1
    switch par.credit_parameter
        case 'CHI_PI' % solve for CHI_PI
            CHI_k   = par.CHI_k;
            CHI_PI  = x(4);
        case 'CHI_k'
            CHI_k   = x(4);
            CHI_PI  = par.CHI_PI;
    end
    CHI_q   = x(5);
    XI      = x(6); 

elseif cali_step == 2
    CHI_k   = par.CHI_k;
    CHI_PI  = par.CHI_PI;
    CHI_q   = par.CHI_q;
    XI      = par.XI;
    c_mu    = x(4); %%%% log normal mean
elseif cali_step == 0
    CHI_k   = par.CHI_k;
    CHI_PI  = par.CHI_PI;
    CHI_q   = par.CHI_q;
    XI      = par.XI;
    c_mu    = par.c_mu;
    c_sig   = par.c_sig;
end
ALPHA    = par.ALPHA; 
CHI_PI_k = CHI_PI;



%----------      Main Equilibrium Conditions   ----------------------------
Z_K    = L * (1 - par.tau_k) * B + (1 - par.DELTA) * (1 - CHI_q - CHI_k); % Z / K

w      = par.ETA / (B / par.A / (1 - par.ETA))^((1 - par.ETA) / par.ETA);

% Functions always have seller' productivity in the 1st argument
% Cutoff value indicates when the division line reaches the upper bound of buyer's productivity

%%%% for capital: the integral
%%%% The range below takes into account the boundaries of productivities
a_b = par.a; % a_b is the persistent level of productivity as a buyer
a_s = par.a; % a_s is the persistent level of productivity as a seller
q_capital = integral2(@(eps, eps_t)fun.capital_gain(eps, eps_t, a_s, a_b),...
fun.l_bound(a_s), min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)),...
@(eps)eps, @(eps)fun.threshold(eps, L, CHI_PI, CHI_PI_k, 1),...
'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

%%%% for money: the integral
a_b = par.a; 
a_s = par.a;
int_l = min(max(L / par.THETA - (1 - par.THETA - CHI_PI - CHI_PI_k) / par.THETA *  fun.h_bound(a_b), fun.l_bound(a_s)), fun.h_bound(a_s)); %need an explanation
int_h = min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)); %need an explanation
q_money   = integral2(@(eps_t, eps)fun.money_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
int_l, int_h, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
+           integral2(@(eps_t,eps)fun.money_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
int_h, fun.h_bound(a_s), @(eps_t)eps_t, @(eps_t)eps_t * 0 + fun.h_bound(a_b),...
'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

%%%% for Holmostrom - Tirol benefit
q_profit   = integral2(@(eps_t, eps) fun.profit_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
int_l, int_h, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
+           integral2(@(eps_t, eps) fun.profit_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI),...
int_h, fun.h_bound(a_s), @(eps_t)eps_t, @(eps_t)eps_t * 0 + fun.h_bound(a_b),...
'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

% Now, assign values to the residules of EE equations

%%%% The residuals of the capital Euler equation. See eq (24)
fval(1) = (1 / par.BETA - (1 - par.DELTA)) / (1 - par.tau_k)...
    - B * (par.q_average_prod + ALPHA_0 * ALPHA * (1 - par.THETA) * q_capital...
    + ALPHA_0 * ALPHA * par.THETA * CHI_k * (1 - par.DELTA) * q_money...
    + ALPHA_0 * ALPHA * par.THETA * CHI_PI_k * q_profit);

%%%% The residuals of the money Euler equation. See eq (25); Below is
%%%% modified version that takes
fval(2) = par.iota - (1 - par.tau_k) * B  * ALPHA * par.THETA * q_money; 



%----------      Other Equilibrium Conditions   ----------------------------
% other eqm variables; only compute if eqm_all = 1

if eqm_all == 1
    c_star  = (XI / (1 - par.tau_h) / w) ^ (1 / (-par.SIGMA));

    q_liq   = integral2(@(eps_t, eps)fun.liq_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
    int_l, int_h, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
    +         integral2(@(eps_t,eps)fun.liq_gain(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K),...
    int_h, fun.h_bound(a_s), @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b),...
    'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

    % K_bar_K = par.q_average_prod + ALPHA_0 * ALPHA * q_capital...
    %     + ALPHA_0 * ALPHA * q_money * (Z_K + (1 - par.DELTA) * CHI_k + CHI_PI * (1 - par.tau_k) * B * par.q_average_prod); 
    
    K_bar_K = par.q_average_prod + ALPHA_0 * ALPHA * q_capital + ALPHA_0 * ALPHA * q_liq; 
     
    a_b = par.a; 
    a_s = par.a;
        
    E_a_K   = ALPHA_0 * ALPHA * integral2(@(eps_t, eps)fun.reall(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, 1),...
            fun.l_bound(a_s), min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)), ...
            @(eps_t)eps_t, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1) ,'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);
    
    E_p_K   = ALPHA_0 * ALPHA *...
            ( integral2(@(eps_t,eps)fun.reall(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, 1),...
            int_l, int_h, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
            + integral2(@(eps_t,eps)fun.reall(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z_K, 1),...
            int_h, fun.h_bound(a_s), @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b), 'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)); 
    
    %%%% note: c + G = Y + par.DELTA * K, and G / Y = par.G_Y
        
    if cali_step ==1
        K      = c_star / (B * K_bar_K / (1 - par.ETA)* (1 - par.G_Y) - par.DELTA);
        Y      = B * K_bar_K / (1 - par.ETA) * K;
        G      = Y * par.G_Y;
    else
        G      = par.G;
        K      = (c_star + par.G) / ( B * K_bar_K / (1 - par.ETA) - par.DELTA);
        Y      = B * K_bar_K  / (1 - par.ETA) * K;
    end

    Z          = Z_K * K;
    H          =  (Y / (par.A * K_bar_K * K)^(1 - par.ETA))^(1 / par.ETA); 
    
    %Cash_ratio = ALPHA_0 * Z / K;
    
    prob_P   = integral2(@(eps_t,eps)fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
            int_l, min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)),...
            @(eps_t)fun.threshold(eps_t,L, CHI_PI, CHI_PI_k, 1), @(eps_t)eps_t - eps_t + fun.h_bound(a_b),'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol)...
              + integral2(@(eps_t,eps)fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
            min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)), fun.h_bound(a_s),...
            @(eps_t)eps_t, @(eps_t)eps_t - eps_t + fun.h_bound(a_b), 'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);        
    
    prob_A   = 2 * integral2(@(eps_t, eps)fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
            fun.l_bound(a_s), min(max(L / (1 - CHI_PI - CHI_PI_k), fun.l_bound(a_s)), fun.h_bound(a_s)),...
            @(eps_t) eps_t, @(eps_t)fun.threshold(eps_t, L, CHI_PI, CHI_PI_k, 1),'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);        
    
    average_2nd_price = 2 * integral2(@(eps_t, eps) ((1 - par.tau_k) * B...
                * ((1 - par.THETA) * eps + par.THETA * eps_t) + 1 - par.DELTA) .* fun.f(eps_t, a_s) .* fun.f(eps, a_b),...
                fun.l_bound(a_s), fun.h_bound(a_s), @(eps_t)eps_t, fun.h_bound(a_b),'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);
    
    %%%% using debt first in transaction (we have tried using cash first)
    debt = integral2(@(eps, eps_t)fun.debt_fun(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_k, Z_K),...
        fun.l_bound(a_b), fun.h_bound(a_b), @(eps)eps, fun.h_bound(a_s), 'method', int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol);

    %%%% coefficient of variation to compute productivity dispersion
    weight     = ALPHA_0 * ALPHA;
    
    moment_1   = (1 - weight) * exp(par.log_eps_mu + par.log_eps_sig^2 / 2) +...
        weight * integral(@(x) x .* fun.g(x, L, par.a), fun.l_bound(par.a), fun.h_bound(par.a));

    moment_2   = (1 - weight) * ((exp(par.log_eps_sig^2) - 1) * exp(2 * par.log_eps_mu + par.log_eps_sig^2) + (exp(par.log_eps_mu + par.log_eps_sig^2 / 2))^2)+ ...
        weight * integral(@(x) x.^2 .* fun.g(x, L, par.a), fun.l_bound(par.a), fun.h_bound(par.a));
    
    prod_std   = sqrt(moment_2 - moment_1^2) / moment_1; 
end

surplus_entry = ALPHA * integral2(@(eps_t, eps)fun.surplus_buy(eps_t, eps, a_s, a_b, B, CHI_q, CHI_PI, CHI_PI_k, CHI_k, Z/K, 1),...
                fun.l_bound(a_s), fun.h_bound(a_s), @(eps_t)eps_t, fun.h_bound(a_b),...
                'method',int_method, 'RelTol', Rel_Tol, 'AbsTol', Abs_Tol) * K / w * XI;
surplus_net = surplus_entry - par.iota * Z / w * XI;  

if cali_step == 1
    fval(3) = E_p_K - par.P_share_target * (E_p_K + E_a_K);
    fval(4) = (E_a_K + E_p_K) - par.R_ratio_target * (E_a_K + E_p_K + par.DELTA);    
    fval(5) = par.Cash_target * Y - ALPHA_0 * Z;
    fval(6) = H / par.H_target - 1;
elseif cali_step == 2   
    c_sig   =  (log(par.surplus_net) - c_mu) / erfinv(2 * par.ALPHA_0 - 1) / sqrt(2); % at the calibrated par.ALPHA_0 and par.surplus_net
    fval(3) = ALPHA_0 - 0.5 * (1 + erf((log(surplus_net) - c_mu) / c_sig / sqrt(2)));
    fval(4) = E_a_K -  ALPHA_0 * prob_A * (par.aqc_spending + par.aqc_spending * par.iota_per_change * par.aqc_elasticity_target);

elseif cali_step == 0
    fval(3) = ALPHA_0 - 0.5 * (1 + erf((log(surplus_net) - c_mu) / c_sig / sqrt(2)));    
end


%----------       Function Output              ----------------------------

fval = fval * 10000; %a simple way to reduce numerical errors
eqm.L = L;
eqm.B = B;
eqm.ALPHA_0 = ALPHA_0;
eqm.c_star          = c_star;
eqm.K               = K;
eqm.Y               = Y;
eqm.G               = G;
eqm.Z               = Z;
eqm.H               = H;
eqm.K_bar_K         = K_bar_K;
eqm.E_a_K           = E_a_K;
eqm.E_p_K           = E_p_K;
eqm.q_capital       = q_capital;
eqm.q_money         = q_money;
eqm.q_profit        = q_profit;
eqm.q_liq           = q_liq;
eqm.w               = w;
eqm.prob_A          = prob_A;
eqm.prob_P          = prob_P;
eqm.average_2nd_price = average_2nd_price;
eqm.prod_std        = prod_std;
eqm.debt            = debt;
eqm.surplus_entry   = surplus_entry;
eqm.surplus_net     = surplus_net;

end

